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CHAPTER  1.  INTRODUCTION 


There  is  much  interest  in  relative  formations  in  low  Earth  orbit  for  many  different  applications  including 
distributed  aperture  systems  and  communications.  While  the  Clohessy-Wiltshire  Hills  (CWH)  equations 
have  been  in  existence  for  sometime,  it  is  more  recently  that  they  have  been  revisited  to  include 
perturbations  in  their  formulation.  Sabol,  et  al1  developed  a  state  transition  matrix  that  included  the  secular 
effects  of  the  oblateness  of  the  Earth  on  relative  trajectories.  In  Refs.  2  and  3,  methods  were  developed  so 
as  to  create  relative  orbits  that  were  J2  invariant.  These  works  also  developed  the  practice  of  using  mean 
orbital  elements  to  design  relative  orbits.  Impulsive  maneuver  schemes  have  been  developed  in  Refs.  3 
and  4  that  operate  on  differential  orbital  elements.  A  part  of  these  approaches  consists  of  a  two  burn 
sequence  at  periapsis  and  apoapsis.  This  limits  the  number  of  correction  maneuvers  to  one  per  orbit.  There 
are  some  similarities  between  Ref.  5  and  this  work,  but  the  physical  variables  used  in  this  report  make  it 
possible  to  formulate  the  nominal  trajectory  in  terms  of  relative  orbital  elements  that  are  more  useful  for 
design  and  analysis. 

There  are  two  major  aims  of  this  report.  First  is  to  develop  a  way  of  designing  relative  orbits  that  allows 
the  designer  to  be  able  to  specify  aspects  of  the  relative  trajectory  while  trying  to  alleviate  the  effects  of  J2 
on  the  orbit.  Second  is  to  implement  a  maneuver  scheme  that  allows  for  multiple  corrections  per  orbit. 
Most  of  the  work  in  this  area  to  date  does  not  allow  for  this  and  problems  arise  when  tolerances  on 
deviation  are  very  low. 
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CHAPTER  2.  RELATIVE  FORMATIONS 


For  this  paper,  we  are  only  concerned  with  the  motion  of  a  deputy  satellite  in  relation  to  a  chief  satellite  on 
a  circular  orbit.  The  position  of  the  deputy  is  in  respect  to  the  chief  in  the  Local  Vertical-Local  Florizontal 
(LVLFI)  frame.  This  frame  of  reference  is  fixed  to  the  chief  satellite.  The  x  direction  is  along  the  zenith 
direction,  y  is  in  the  alongtrack  direction,  and  z  is  perpendicular  to  the  orbital  plane. 


Using  the  parameterization  from  Ref.  6,  the  motion  of  the  satellite  is  characterized  according  to  the 
following  equations: 


X  =  — ^  COS  P  +  xd 

y  =  ae  sin  /3  +  yd  nxdt 

z  =  zmaxSin  {y  +  P) 


■  ae  ■  a 

x  —  — —  n  sin  p 
2 

„  3 

v  =  aen  cos  p  -  —  nxd 
z  =  zm^n  cos  (y  +  /3) 


(l) 


The  equations  above  describe  a  relative  ellipse  that  travels  along  the  alongtrack  direction.  The  relative 
orbital  elements  are  as  follows:  ae  is  the  major  axis  of  the  ellipse,  xd  and  yd  describe  offsets  to  the  relative 
ellipse,  in  the  zenith  and  alongtrack  directions  respectively,  zmax  is  the  maximum  excursion  out  of  the 
orbital  plane,  y  describes  the  orientation  of  the  ellipse  and 

y 3  =  nt  (2) 

These  parameters  are  depicted  in  Fig.  1  for  xd  =  yd  =0.  The  parameterized  values  for  [)  and  y  are  related  to 


the  geometrical  values  in  the  figure,  /?  and  y  ,  by: 

tan  (3  =  2  tan  (i 
tan  y  =  2  tan  y 

The  projection  of  the  path  of  the  deputy  satellite  on  to  the  x-y  plane  forms  a  2x1  ellipse. 


(3) 


2 


Figure  1:  Relative  Orbit 


For  unperturbed  motion,  ae,  xd,  yd,  zmax,  and  y  are  all  constant.  The  only  time  varying  element  is  /?.  This 

gives  a  set  of  six  elements  that  characterize  a  relative  formation  which  will  be  referred  to  as  the  “relative 
orbital  elements.”  For  a  closed  formation,  i.e.  one  in  which  all  the  satellites  would  remain  in  close 
proximity  rather  than  drift  apart,  the  secular  term  in  the  alongtrack  equations  must  be  eliminated  by  setting 
x,j=0.  Setting  is  equivalent  to  constraining  the  orbit  of  the  chief  and  the  deputy  to  have  the  same  semi 
major  axis.  In  an  unperturbed  state  this  is  also  equivalent  to  setting  the  periods  of  the  two  satellites  equal  to 
each  other.  This  is  the  only  way  to  have  a  closed  formation  for  unperturbed  motion. 
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CHAPTER  3.  OBLATENESS  EFFECTS  ON  ORBITAL  ELEMENTS 


A.  Definition  of  Orbital  Elements 

As  defined  in  Ref.  7,  there  are  six  classical  orbital  elements  which  characterize  the  orbit.  The  first  is  the 
semi  major  axis,  a,  which  relates  the  size  of  the  orbit.  Second  is  eccentricity,  e,  which  determines  the  shape 
of  the  orbit,  i.e.  circular,  elliptical,  parabolic,  or  hyperbolic.  For  the  purposes  of  this  paper,  we  are  only 
concerned  with  circular  orbits  which  have  an  eccentricity  of  zero  and  elliptical  orbits  whose  eccentricities 
are  greater  than  zero  and  less  than  one.  Third  is  the  inclination,  i,  which  relates  the  tilt  of  the  orbit  in 
relation  to  the  equatorial  plane.  Fourth  is  the  Right  Ascension  of  the  Ascending  Node,  Q  ,  which  measures 
the  angle  between  the  line  of  Aries  in  the  inertial  frame  and  the  point  in  the  orbit  where  the  satellite  crosses 
from  the  south  to  the  north  hemisphere,  also  known  as  the  ascending  node.  Fifth  is  the  argument  of 
perigee,  at,  which  measures  the  angle  from  the  node  to  the  point  of  closest  approach  to  the  Earth,  or 
perigee.  Last  is  the  true  anomaly,  f,  which  measures  the  angle  between  perigee  and  the  position  of  the 
satellite. 

For  this  paper,  several  of  these  orbital  elements  are  changed  so  as  to  avoid  singularities  and  other 
complications.  In  a  circular  orbit,  there  is  no  perigee,  so  the  angle  used  to  determine  the  satellites  current 
location  is  the  true  argument  of  latitude,  0,  which  is  defined  as  9=a>+f.  Also  in  a  circular  orbit,  the 
eccentricity  is  zero,  which  can  cause  singularities  in  several  of  the  conversions  to  be  used  later  in  the  paper. 
To  overcome  this,  two  new  variables  are  introduced:  qi  and  q2.  They  are  defined  as6 

qi=e  cosro  and  q2=e  sinco  (4) 

To  determine  the  orbital  elements  of  the  deputy,  the  relative  orbital  elements  are  used  with  Eqs  (1)  to 
determine  the  relative  state  of  the  deputy  in  the  LVLH  frame.  The  chief  s  position  and  velocity  in  the 
inertial  frame  of  reference  are  then  transformed  into  a  rotating  frame  of  reference  using  Euler  angle 
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rotations  as  laid  out  in  Vallado7.  The  position  of  the  deputy  is  found  by  adding  the  relative  position  vector 
to  the  rotating  position  vector  of  the  chief.  The  velocity  of  the  deputy  in  the  rotating  frame  is  found  by 
adding  the  relative  velocity  and  the  velocity  of  the  LVLH  frame  to  that  of  the  chief.  The  position  and 
velocity  of  the  deputy  in  the  rotating  frame  are  then  transformed  into  the  inertial  frame  through  Euler  angle 
rotations.  The  deputy  orbital  elements  are  then  determined  by  using  the  deputy  position  and  velocity 
following  processes  laid  out  in  Ref.  7. 


B.  Equations  of  Motion 

To  numerically  propagate  an  orbit,  the  equations  of  motion  need  to  be  developed.  From  Ref.  7,  the  basic 
two-body  equation  of  motion  is 


Separating  this  equation  into  its  components  gives 


(5) 


x  =  — Tx  y  = 

r  r 


H  ..  H 

— y  z=--rz 


(6) 


This  formulation  is  the  unperturbed  two-body  equations  of  motion.  To  include  the  perturbations  due  to  J2, 
we  take  the  partial  derivatives  of  the  disturbing  function  and  add  these  effects  to  Eqs  (6).  Doing  this  yields 
the  following  equations  of  motion  used  to  numerically  propagate  the  satellite’s  orbit. 


u  3  J-,R2,x^ 

X  =  — ~rX - ~  / 

r3  2  r5 


1- 


5z 


2  T 


M  3  J2Rey 

y=-—y — H- 

r  2  r 


1- 


5  z 


2  \ 


V 


(7) 


H  iJ2R2.^ 

Z  = - rZ - - - 

r  2  r~ 


3- 


5z 


2  3 


V 


C.  Oblateness  Effects 

When  considering  unperturbed  two  body  motion,  the  only  time  varying  orbital  element  is  true  anomaly. 
This  is  not  the  case  when  perturbations,  especially  the  oblateness  of  the  Earth,  are  included  in  the  equations 
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of  motion.  When  the  oblateness  of  the  Earth,  specified  by  the  J2  term  in  the  harmonic  expansion  in  Ref.  7, 
is  included,  all  of  the  orbital  elements  experience  some  periodic  oscillations  in  their  values.  This  is  most 
easily  seen  in  the  values  for  semi  major  axis  and  inclination.  An  example  case  is  illustrated  using  the 
following  initial  orbital  elements  for  a  satellite: 

a=  7378  km,  0=0,  i=50  deg,  qi=0,  q2=0,  Q=0  deg 


The  orbit  of  the  satellite  is  propagated  using  Eqs  (7)  and  a  numerical  integration  in  MATLAB.  Time 
histories  of  the  semi  major  axis  and  inclination  are  shown  in  Figures  2  and  3,  respectively,  for  a  duration  of 
10  orbits.  From  these  graphs,  it  is  apparent  that  the  mean  values  are  not  the  initial  values.  The  semi  major 
axis  average  value  is  just  under  7373  kilometers,  from  observation.  Because  of  this  difference  between  the 
initial  values  and  the  average  value,  we  will  define  two  different  sets  of  orbital  elements.  The  osculating 
orbital  elements  will  be  the  instantaneous  orbital  elements  at  any  point  in  time.  The  mean  orbital  elements 
will  be  the  average  orbital  elements  over  time. 
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Figure  2:  Time  History  of  Semi  Major  Axis 


6 


50.005 


50 


49.995 


49.99 


.2  49.985 


^  49.98 


49.975 


49.97 


49.965 L 


0 


l_I 


4  5  6 

time  (orbits) 


10 


Figure  3:  Time  History  of  Inclination 


D.  Characterizing  Perturbation  Effects 

To  set  up  initial  conditions  that  result  in  desired  mean  elements  a  conversion  is  used.  To  develop  this 
conversion,  the  potential  for  J2  gravity  perturbation  needs  to  be  reformulated  to  show  its  effects  on  orbital 
elements.  Following  Kozai8,  the  potential  is 

R(r)  =  ^1fe  (_  j.sin2  ^  +  1)  (8) 

r  2  2 

where  (f)  is  the  latitude  of  the  sub-satellite  point.  Using  relationships  from  spherical  trigonometry,  we  are 
able  to  rearrange  the  equation  into  the  form 

R  =  ^^-[l--sm2i(l-cos20)]  (9) 

2  r  2 
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Kozai  then  breaks  this  disturbing  potential  into  four  parts:  first  and  second  order  secular  terms  (linear  in 
time),  long  period,  and  short  period.  To  break  out  the  first  order  secular  effects,  Rb  a  Fourier  series  is  used 
to  isolate  those  terms.  A  relationship  from  Tisserand9  is  required  and  yields  first  order  secular  effects  of 


*.= 


rt 

2  a3  1 


^-sin2  i)(l-e2)  ^ 


(10) 


As  an  example,  to  find  the  first  order  secular  effects  on  RAAN,  this  equation  is  inserted  into  Lagrange’s 
planetary  equations  to  find 


n=- 


1  SR,  _  3 fjJ2R) 


nabsini  Si  2 p1ain 


cos  i 


(11) 


where  n  is  the  perturbed  mean  motion. 


We  also  need  to  determine  the  short  period  effects  on  orbital  elements.  Using  the  potential  equations,  the 
short  period  effects  are 

R4=R-R,  (12) 

Again,  as  an  example,  we  will  use  the  effect  on  RAAN.  The  disturbing  potential  is  placed  into  the  equation 


n  =  1 

2  •  •  C*. 

na  sin  i  oi 


(13) 


only  retaining  the  terms  to  0{e)  .  This  equation  needs  to  be  integrated  to  determine  the  effects.  First  we 
need  to  change  the  variable  of  integration  from  time  to  true  anomaly  using  the  relationship 


St  =  - 
n 


r 

\a  J 


1 


a/i  —  < 


df 


(14) 


Using  this  equation  we  are  able  to  rearrange  the  terms  to  get 


r  3J0i?2cosz 

f  a 'l 

f  r'] 

2 

2  e 

(1-COS261) 

~ 

J  2  a2 

Kr  J 

l a) 

8f 


(15) 
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Before  we  integrate  we  need  to  get  this  integral  equation  into  terms  of  true  anomaly.  Using  a  Taylor  series 
expansion  and  the  formula 


r  = 


a(  l-e2) 


(16) 


1  +  ecos  / 

along  with  several  trigonometric  identities  we  are  able  to  get  Eq  (15)  into  a  form  suitable  for  integration. 
Integrating  and  rearranging  the  terms  leads  to  the  form 


^  3J7R~  cosi 

Q  = - - — £ - 

2a2 


3esin /-  —  sin2<9-  —  esin(/  +  2co)-  —  esin(3/  +  2a>) 
2  2  6 


(17) 


The  same  process  is  followed  for  the  other  orbital  elements. 

To  transform  a  set  of  elements  from  mean  to  osculating  the  following  method  is  used: 

easc  =  emean  +  s(Ae,p  +  Aespl  +  Aesp2)  (18) 

Where  e  is  \aOi  qt  q2  and  £  =  — J2R 2 .  The  long  period  variation  due  to  Earth  oblateness,  Aelp  , 

and  the  short  period  variations,  Aesp'  ,Aespl ,  can  be  determined  from  the  preceding  development.  The 
exact  formulation  used  for  this  paper  is  from  Alfriend  and  Gim10,  based  on  work  by  Brouwer11  and 
Lydanne12.  The  equations  used  are  listed  in  Appendix  I. 


E.  Secular  Effects 

In  addition  to  the  periodic  perturbations,  the  Mean  anomaly,  M,  the  Argument  of  Perigee,  ro,  and  the 
RAAN  experience  secular  changes.  These  angular  rates  are  obtained  from  general  perturbation  theory: 


(n  \2 

Q  =  -1.5  J,n\ 


cost 


V  P  J 

R 


(  d  T2 


d)  =  0.75  J7n 

PJ 

M  =  n  +  0.75  J2n\ 


(5  cos2  /  -1 ) 


(19) 


fR  ^ 


\  P  J 


V 1  -  e 1  (3 cos2  i  —  1 
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Therefore  the  time-varying  mean  elements  are 


a(t )  =  a0 
0(t)  =  0Q  +  6t 
i(t)  =  i0 

(20) 

?2(0  =  qi„ 

O(f)  = 

As  described  above,  the  osculating  elements  are  obtained  by  adding  the  periodic  terms  to  these  mean 
elements. 

For  this  paper,  an  orbital  period  is  defined  as  the  time  it  takes  for  the  deputy  to  travel  through  2  n  radians 
in  the  precessing  orbital  plane.  As  such,  we  will  also  define  the  mean  motion  of  the  satellite  to  be  its 
perturbed  mean  motion,  in  other  words,  the  mean  motion  will  be  equal  to  the  unperturbed  mean  motion 
plus  the  secular  rate  of  change  of  the  mean  anomaly  due  to  J2,  as  seen  in  Eqs  (19).  This  allows  for  the 
mean  motion  to  account  for  the  secular  changes  to  the  Mean  anomaly  due  to  Earth  oblateness. 
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CHAPTER  4.  OBLATENESS  EFFECTS  ON  RELATIVE  TRAJECTORY 


A.  Effects  on  an  In-plane  Formation 

The  largest  effect  of  the  oblateness  perturbation  is  an  induced  drift  in  the  relative  trajectory.  To  see  this,  a 
relative  formation  is  set  up  with  one  chief  and  one  deputy.  The  initial  orbital  elements  of  the  chief  are: 

a=  7378  km,  9=0,  i=50  deg,  qi=0,  q2=0,  0=0  deg 
The  initial  relative  orbital  elements  for  the  deputy  are: 

ae=  500  m,  xd=0,  yd=0,  zmax=0,  y=0  deg,  (1=0  deg 

This  is  considered  to  be  an  “in-plane”  formation  because  the  relative  ellipse  is  in  the  plane  of  the  orbit. 
This  is  indicated  by  zmax=0.  These  conditions  are  used  to  determine  the  initial  conditions  of  the  deputy 
using  Eqs  (1)  and  the  processes  outlined  in  Chapter  2.  The  states  of  the  two  satellites  are  then  propagated 
forward  in  time  using  Eqs  (7),  which  include  oblateness  perturbations.  A  Runge-Kutta  variable  step  size 
numerical  integrator  in  MATLAB  is  used  for  the  simulation.  The  results  are  shown  in  Figure  4  for  a  time 
period  of  10  orbits  and  indicate  a  drift  on  the  order  of  five  meters  per  orbit.  This  is  problematic  for  a 
station-keeping  scheme  in  that  the  satellite  would  constantly  be  fighting  drift  to  remain  in  a  nominal 
relative  orbit. 
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Figure  4:  Relative  Formation  Drift 


B.  Effects  on  an  Out-of-plane  Formation 

This  effect  is  also  seen  in  an  out-of-plane  formation.  To  illustrate  this,  the  same  initial  conditions  for  the 
chief  will  be  used.  For  the  deputy,  all  relative  orbital  elements  will  remain  the  same  with  the  exception  of 
Zmax*  For  this  example,  zmax=500  m  to  achieve  out-of-plane  relative  motion.  The  results  are  seen  in 
Figure  5.  This  formation  has  a  drift  on  the  order  of  fourteen  meters  per  orbit. 
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y  (m)  X  (m) 

Figure  5:  Oblateness  Effects  on  Out-of-Plane  Motion 
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CHAPTER  5.  J2  INVARIANT  RELATIVE  ORBITS 


A.  Orbital  Element  Initialization 

Eliminating  drift  in  relative  orbits  requires  several  steps.  The  first  step  is  to  consider  the  initial  orbital 
elements  of  both  the  chief  and  the  deputy  to  be  mean  elements  and  then  convert  them  to  osculating 
elements.  Doing  this  allows  for  the  periodic  oscillations  of  the  elements  to  be  taken  into  account  when 
initializing  a  formation.  Care  needs  to  be  taken  in  choosing  which  method  to  use  to  transform  your 
elements.  Small  errors  can  lead  to  bad  initial  conditions  for  the  formation.  The  process  used  for  this  paper 
is  outlined  in  Chapter  II.  For  the  rest  of  this  paper,  when  discussing  formation  design,  all  orbital  elements 
are  assumed  to  be  mean  orbital  elements. 


B.  Matching  Periods  and  Nodal  Precession  Rates 

As  mentioned  in  Chapter  III,  Earth  oblateness  causes  secular  change  in  an  orbit’s  mean  anomaly,  argument 
of  perigee,  and  RAAN.  Differences  in  these  rates  between  a  chief  and  a  deputy  cause  unwanted  drift  in  a 
formation.  To  negate  this  effect,  we  perturb  the  deputy’s  orbit  to  match  the  period  and  nodal  precession 
rate.  These  two  conditions  are  what  are  referred  to  as  the  J2  invariant  orbit  in  Ref.  2.  They  are  expressed 
as: 


Ci 


c 


(21) 


0C  -  Mc  +  d>c  -0d 

The  subscripts  of  “c”  and  “d”  refer  to  the  chief  and  deputy,  respectively.  Formulas  for  the  angular  rates  are 
from  Eq  (18). 


To  match  the  period  and  nodal  precession,  we  take  the  variation  of  Eqs  (19).  Problems  arise  for  circular 
chief  orbits  because  the  terms  accompanying  eccentricity  will  drop  out  of  the  equations.  Therefore,  as  in 
Ref.  13,  we  define  the  new  variable: 
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(22) 


7  = 


Also,  for  convenience,  we  define 


3J2Re2n 

1  2  4 

2  a  r) 


(23) 


Taking  the  variation  of  Eqs  (19)  gives  us: 

8{M  +  oj)  =  s[n  +  C(1  - 1 sin 2  i)rj  +  C(fcos2  i  -|)J 

=  -  £  [3/i  +  7C(1  - 1  sin2  07  +  7C(|  cos2  i  - 1)]& 

-  c[3(l  -  4 sin2  0  +  f  (f  cos2  i  -  4)]fS7  -  C  cos  i  sin  i(3rj  +  5 )8i  ^ 

=  S(-  C  cos  i)  =  cos  iSa  +  cos  idr/  +  C  sin  idi 

where  the  orbital  elements  are  mean  chief  elements.  Using  Eqs  (24),  it  is  then  possible  to  determine  a 
relative  deputy  orbit  that  will  not  drift. 


C.  Designing  a  Relative  Orbit 

To  aid  in  designing  a  relative  trajectory  that  does  not  drift,  we  have  reformulated  Eqs  (24)  in  terms  of 
relative  orbital  elements.  To  do  this,  we  first  need  to  take  the  variation  of  Eq  (22).  We  define 


§T]  =  rjd  —  rjc  (25) 

Then  we  use  Eq  (22)  and  the  fact  that  our  chief  is  circular  to  get 

Srj  =  ^l-el  -1  (26) 

Since  our  chief  is  circular  we  can  use  the  substitution 

ed=8e  (27) 


to  change  Eq  (26)  to  the  form 

Sri  =  ^  1-&2  - 1  (28) 

and  finally  use  a  binomial  series  expansion  to  reach  the  equation 

S?l~-^Se2  (29) 

We  then  make  use  of  a  linear  conversion  between  orbit  element  differences  and  relative  orbit  elements  of 
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da  =  xd 

Si  =  cos(d  -(/  +  /?)) 
a 

Se  = ^ 

2  a 

to  get  Eqs  (24)  in  terms  of  relative  orbit  elements. 

8  (M  +  (b)  =  -  \$n  +  7  C(1  - f  sin2  i)rj  +  7C(4  cos2  z  - 

+  N1  ~  f  sin2  0  +  7  (I COs2  *  “  1)]  ae2 

-  £  cos  z  sin  i  (3/7  +  5)  zmax  cos[<9 -(/  +  ^)] 

Stl^^cosiXj-^-r cos /  a,2  +  f  sin z  zmax  cos[# -  (7  +  /?)] 

This  formulation  allows  for  one  of  the  set  of  relative  elements  and  element  combinations  {xd,  ae2, 
zmaxCOs[9-(y+P)] }  to  be  fixed  and  determine  the  other  two  from  the  linear  system  of  equations  above. 
These  T  invariant  relative  orbit  elements  can  then  be  used  to  determine  the  orbital  elements  of  the  deputy. 

For  example,  the  designer  may  wish  to  consider  various  values  of  y  and  look  at  the  values  of  xd  and  ae  take 
on  for  J2  invariance.  These  relationships  are  shown  in  Figure  6  for  a  zmax  of  50  meters.  Note  that  for  a 
value  of  n  1 2  that  ae  becomes  very  small.  This  is  because  there  is  no  inclination  difference  between  the 
deputy  and  the  chief  for  this  value  of  y . 


(30) 
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Figure  6:  Conditions  for  J2  Invariance  with  Varying  y 

For  an  out-of-plane  formation,  problems  occur  due  to  the  semi  major  axis  squared  in  the  denominator  of  the 
ae  term  in  Eqs  (31).  This  makes  the  coefficient  of  ae  very  small,  so  large  changes  in  ae  are  necessary  of  any 
non-negligible  change  in  the  other  relative  orbital  elements.  Thus,  matching  both  the  period  and  the  nodal 
rates  for  significant  zmax  with  relatively  small  ae  values  (kilometers  or  less)  is  physically  impossible.  As 
seen  if  Figure  6,  for  relatively  small  values  of  zmax,  the  value  of  ac  increases  dramatically.  Therefore,  in  this 
report  we  will  only  match  the  period  between  the  deputy  and  the  chief.  The  repercussions  of  this  will  be 
seen  in  the  following  section. 
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Figure  7:  Various  Conditions  for  J2  Invariance  with  Varying  zmax 

To  demonstrate  the  effectiveness  of  period  and  nodal  rate  matching  we  consider  both  an  in-plane  and  out- 
of-plane  example.  For  the  in-plane  example,  solving  Eqs  (31),  gives  the  following  relative  orbital  elements 
for  the  deputy: 

ae=  500  m,  xd=-3.7  x  10"5  m,  yd=0,  zmax=3  x  10"16,  y=0  deg,  P=0  deg 
Notice  that  the  values  for  xd  and  zmax  are  very  close  to  the  original  values.  Because  there  is  no  inclination 
difference  between  the  chief  and  the  deputy,  the  period  and  nodal  rate  matching  conditions  are  very  nearly 
met  by  simply  setting  xd=0. 

These  equations  are  use  till  in  designing  and  analyzing  a  relative  orbit.  However,  using  these  linear 
conversions  to  initialize  an  orbit  leads  to  unnecessary  error  in  the  solution.  Instead,  Eqs  (21)  are  used 
directly  to  initialize  the  orbit.  Doing  so  leads  to  the  following  initial  orbital  elements  for  the  deputy: 

a=  7378999.999963  m,  0=0  deg,  i=50  deg,  qi=3.389  x  10'5,  q2=0,  Q=0  deg 
For  these  orbital  elements  to  produce  a  J2  invariant  orbit,  they  are  considered  to  be  mean  orbital  elements. 
Utilizing  the  mean  to  osculating  conversion  gives  the  following  initial  osculating  orbital  elements  for  the 
chief  and  the  deputy: 

ac=  7383251.179  m,  0C=  0  deg,  ic=  50.0171  deg,  qic=  7.384  x  10"4,  q2c=0  ,  Qc=0  deg 
ad=  7383251.786  m,  0d=  0  deg,  id=  50.0171  deg,  q[d=  7.723  x  10"4,  q2d=0  ,  Qd=0  deg 
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Using  the  MATLAB  simulation  previously  described,  the  relative  trajectory  created  by  using  these  initial 
conditions  is  seen  in  Figure  8.  The  drift  from  Figure  4  is  completely  negated  by  the  mean  to  osculating 
conversion.  Over  the  100  orbits  graphed,  the  drift  in  the  formation  is  only  .9  m,  or  9  millimeters  per  orbit. 


Figure  8:  In-plane  J2  Invariant  Relative  Trajectory  (100  Orbits) 

For  the  out-of-plane  example,  as  mention  above  we  cannot  find  a  fully  J2  invariant  orbit,  so  we  only  match 
the  period.  Specifying  two  of  the  three  relative  orbit  parameters  in  the  period  matching  equation  in  Eqs 
(31)  will  allow  for  an  out-of-plane  formation  that  has  reduced  drift  versus  its  non-matched  counterpart. 
Using  this  method  for  an  out-of-plane  formation  yields  the  following  initial  relative  orbital  elements: 

ae=  500  m,  xd=  -1.590  m,  yd=0,  zmax=500  m,  y=0  deg,  P=0  deg 
The  perturbation  to  xd  is  much  larger  than  that  of  the  in-plane  formation.  Using  these  initial  conditions  to 
simulate  the  relative  formation  leads  to  the  relative  trajectory  in  Figure  9.  The  drift  of  this  formation  is 
greatly  reduced  from  that  of  the  trajectory  in  Figure  5. 
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Figure  9:  Out-of-Plane  Period  Matching  -  10  Orbits 


This  reduction  in  drift  from  an  unmatched  trajectory  will  keep  the  formation  from  fighting  the  full  effects 
of  drift  caused  by  J2.  Even  so,  there  is  still  drift  in  the  formation  of  about  1.84  meters  per  orbit.  This  is 
caused  by  the  unmatched  differential  secular  rate  of  change  between  the  RAAN  of  the  chief  and  the  deputy. 
The  differential  nodal  rate  causes  in-plane  drift  to  occur  as  discussed  in  the  next  section. 


D.  Drift  Due  to  Differential  Nodal  Rates 

To  characterize  this  drift,  assume  that  the  chief  and  the  deputy  are  collocated.  We  can  do  this  without  loss 
of  generality.  With  a  nearly  circular  chief  and  period  matching  conditions  enforced,  the  chief  and  the 
deputy  will  both  recross  the  equatorial  plane  at  the  same  time  one  period  later.  Over  that  period  of  time,  T, 
a  difference  in  the  RAAN  of  the  two  satellites  will  have  occurred.  The  distance  between  the  chief  and  the 
deputy  will  be  given  by 

(Qd-Qc)Tac  (32) 

The  component  parallel  to  the  orbit  plane  (see  Figure  10)  is: 

Drift  =  AClTa  COS  i  (33) 
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Drift  =  AQ  T  a  cos  i 


Figure  10:  Drift  Due  to  Differential  Nodal  Rate 

Using  the  example  from  Figure  9,  the  calculated  drift  woidd  be 

In-plane  drift  =  (Q^  —Clc)Tac  COS ic  =  1.88  m 

which  matches  up  nicely  with  the  observed  drift  seen  in  Figure  1 1  of 

In-plane  drift=(517.3-500.7)/9=1.84  m 


Figure  11:  Blow  up  of  Figure  9  to  Determine  Drift 


(34) 


(35) 


In  addition  to  the  in-plane  drift,  there  is  also  some  out-of-plane  motion  associated  with  nodal  drift.  This 
motion  affects  y  and  zmax. 


In  the  J2  invariant  work  by  Schaub  and  Alfriend2,  their  process  for  finding  an  invariant  relative  orbit 
requires  defining  one  of  three  deputy  orbital  elements,  either  semi  major  axis,  eccentricity,  or  inclination. 
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and  then  finding  the  other  two  to  initialize  the  relative  orbit.  The  problem  with  this  approach  is  that  the 
relative  orbital  elements  that  define  such  an  orbit  are  unable  to  be  changed  and  do  not  always  form  a 
relative  orbit  that  is  desired. 

For  formation  design  purposes,  this  approach  lacks  the  ability  to  define  the  J2  invariant  orbit  that  meets  the 
requirements  of  the  formation  design.  The  approach  laid  out  here  allows  for  a  relative  orbit  that  not  only 
has  close  to  the  desired  relative  orbital  elements,  but  also  reduces  drift  to  an  amount  that  can  be  corrected 
by  a  stationkeeping  scheme. 
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CHAPTER  6.  GUIDANCE  APPROACH 


A.  Defining  the  Nominal  Trajectory 

For  a  stationkeeping  scheme,  two  different  nominal  trajectories  need  to  be  defined:  the  trajectory  used  to 
monitor  deviation  and  the  target  state  for  a  maneuver.  Deviations  from  the  nominal  will  be  measured  in  the 
relative  state.  The  least  computationally  intensive  solution  is  to  use  Eqs  (1).  This  is  the  trajectory  we  will 
use  to  monitor  deviation.  When  we  initialize  the  formation,  we  perturb  the  deputy  to  create  a  period 
matched  orbit.  When  using  the  parameters  to  monitor  the  formation,  we  will  set  xd=0,  since  the  drift  has 
ideally  been  removed.  The  mean  motion  used  will  be  that  of  the  J2  perturbed  mean  motion.  To  propagate 
forward  our  relative  trajectory,  we  will  propagate  forward  /?to  a  time  T  by  using: 

P{T)  =  P0+MT  (36) 

For  the  target  state  at  the  end  of  a  maneuver,  we  will  use  a  slightly  different  approach.  For  target 
conditions,  we  need  to  take  into  account  targeting  mean  elements  (versus  osculating  elements)  and  we  will 
need  to  perturb  the  target  to  remove  any  unnecessary  drift.  We  will  first  propagate  forward  the  mean 
argument  of  latitude  of  the  chief.  Likewise,  the  chief  s  anomaly  is  propagated  forward  by  simply  using  the 
perturbed  mean  motion: 

0(T)  =  0o+MT  (37) 

In  the  relative  state,  the  RAAN  is  ignorable,  so  it  is  not  updated.  The  chief  orbital  elements  and  the  relative 
orbital  elements  of  the  deputy  at  time  T  are  then  used  to  find  the  desired  orbital  elements  of  the  deputy. 

To  remove  unwanted  drift,  the  deputy’s  mean  orbital  elements  are  then  perturbed  using  the  period  matching 
process  described  in  Chapter  VI.  The  desired  orbital  elements  of  the  chief  and  deputy  are  then  converted  to 
osculating  elements.  Once  these  are  found,  the  target  state  of  the  deputy  is  found  in  a  LVLH  frame  in  a 
process  similar  to  that  mentioned  in  Chapter  II,  found  in  standard  astrodynamical  texts7.  It  is  this  relative 
position  and  velocity  that  are  used  for  targeting  the  maneuver. 
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B.  Calculating  AV  for  Maneuvering 


The  nominal  trajectory  used  to  determine  when  the  deputy  needs  to  maneuver  is  the  nominal  position  given 
by  Eqs  (1).  To  measure  deviation,  the  nominal  relative  position  is  subtracted  from  the  actual  relative 
position  and  the  magnitude  of  this  position  vector  is  considered  the  deviation.  A  radius  of  a  sphere,  or 
deadband,  is  specified  that  we  want  the  satellite  to  remain  inside.  To  determine  if  a  maneuver  is  warranted, 
we  check  the  deviation  to  see  if  it  has  reached  a  certain  percentage  of  that  deadband.  If  it  is  at  that  distance 
or  further,  we  enact  a  maneuver.  Each  stationkeeping  maneuver  consists  of  two  burns.  This  gives  six 
degrees  of  freedom  to  completely  match  the  nominal  state.  The  first  maneuver  will  send  the  satellite  back 
to  a  nominal  position,  &(T)  ,  where  T  is  the  time  at  the  end  of  the  maneuver.  The  maneuver  needs  to  be 
short  enough  that  several  maneuvers  can  take  place  per  orbit. 

To  determine  the  trajectory  that  will  return  our  satellite  to  nominal,  we  will  need  to  develop  a  state 
transition  matrix  for  the  LVLH  frame  of  reference.  We  follow  the  process  laid  out  in  Pmssing  and 
Conway14.  Developing  the  necessary  equations  requires  us  to  start  in  an  inertial  frame  of  reference,  seen  in 
Figure  12. 


Figure  12:  Inertial  Frame  of  Reference 
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Using  the  equation  of  motion  in  a  general  gravity  field, 


r  =  g(r ) 


(38) 


along  with  the  relationship  of 


r=r*+SrI  (39) 

we  are  able  to  get  a  series  expansion  of  the  equation  of  motion  in  a  general  gravity  field 

s\ct( 

r*+SrI=g(r*+SrI)  =  g(r*)-\ - -Srj+0  (40) 

Sr* 

where  the  higher  order  terms  are  ignorable  due  to  the  large  difference  in  size  between  the  position  of  the 
chief  and  SFj .  Using  this  equation  and  the  fact  that  the  reference  orbit  satisfies  the  equations  of  motion 
gives  us  a  linear  differential  equation  for  our  relative  trajectory  given  in  Eq  (41). 

Sr,  =  &I  =  G{r*)SrI  (41) 


To  make  this  equation  specific  to  our  application,  we  use  the  inverse  square  gravitational  force  equation 

G(r)  =  ^j (3 rrT  - r2I3 )  (42) 

r 

along  with  the  relationships 


2  H 
,  n  = 


*3 


(43) 


to  get  the  specific  equations  of  motion  in  an  inertial  frame  of 


Srj  =  n‘ 


2  0  0 
0-10 
0  0-1 


Sr, 


(44) 


Next,  we  need  to  convert  this  equation  into  the  rotating  LVLH  frame.  To  do  this,  we  will  use  the 
acceleration  expansion 

Sr  =  SrI  -IcoxSr  -cox-Stj  -  x  (<y  x  (3F, )  (45) 

to  get  the  equations  of  motion  in  the  proper  frame.  The  angular  rate,  CO  ,  is  defined  as 
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CD  =  \0 


Assuming  constant  mean  motion  gives 

x  =  3  n2x  +  2  ny 

y  =  -2nx 

2 

z  =  —n  z 

These  equations  are  the  Clohessy- Wiltshire  equations.  To  get  the  state  transition  matrix  for  relative 
position  and  velocity,  these  differential  equations  are  integrated  to  get 


, ,  , .  _  .  sin  nt  .  2 .  . 

x{t)  =  (4-3cos«T)*o H - xo  +— (1-coshT)v0 

n  n 


4  sin  nt-3nt 


y(t )  =  6(sin  nt  -  nt)x0  +  y0  —  (1  -  cos  nt)x0  H - 

n  n 

z(t)  =  z0  cos  nt  +  —  sin  nt 
n 

x(t)  =  (3/7  sin  nt)x0  +  (cos  nt)x  0  +  (2  sin  nt)y0 
y(t)  =  -6(1  -  cos  nt)x 0  -  (2  sin  nt)x 0  +  (4  cos  nt  -  3) y0 
z(t)  =  -z0n  sin  nt  +  z0  cos  nt 


The  terms  are  then  put  into  a  state  transition  matrix  with  the  form 


X(t)  =  ®X(t0) 


where 


This  state  transition  matrix  can  then  be  partitioned  into  the  form 


0  0 

rr  rv 

0.  .  0  , 


These  partitions  can  then  be  used  to  determine  what  velocity  would  be  needed  to  maneuver  to  a  certain 
position  after  a  certain  time  period  by  the  equation 

AV,  =  ©^(r)-1^-©^^)]-^)  (52) 
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where  the  ®  matrices  are  partitions  of  the  state  transition  matrix  and  c>r(0)  and  c)v(O)  are  the  relative 


position  and  velocity  at  the  beginning  of  the  maneuver.  The  position  at  the  end  of  maneuver,  dF(T )  ,  is 
determined  using  the  process  laid  out  at  the  beginning  of  this  chapter. 

The  second  burn  in  the  stationkeeping  maneuver  removes  any  velocity  relative  to  the  reference  so  that  the 
deputy  stays  on  the  nominal  trajectory  is  simply 

A  V2=Svnom(T)-Sv(T)  (53) 

where  ( T)  is  determined  from  the  process  at  the  beginning  of  this  chapter. 

Due  to  the  assumption  that  the  satellite  will  have  three  axis  capability  and  will  not  have  to  slew  to  perform 
a  maneuver,  the  total  A  V  is  given  by 

A(,  =  ZZ|AFi/j  <*> 

*=  1  7=1 

where  i  refers  to  the  first  and  second  burns  and  j  refers  to  the  X  ,  y  ,  and  2  direction. 
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CHAPTER  7.  MONTE  CARLO  ANALYSIS 


A.  Monte  Carlo  Setup 

To  evaluate  the  performance  of  the  stationkeeping  scheme,  a  Monte  Carlo  analysis  is  performed  for  100 
runs  of  20  orbits  (2000  total  orbits).  The  chief  satellite  is  placed  on  a  perfectly  known  orbit  with  a  mean 
eccentricity  of  zero.  The  deputy  satellite  is  initially  placed  on  the  desired  relative  formation  but  then 
corrupted  using  the  navigation  uncertainties.  The  chief  and  deputy  satellite  dynamics  are  then  simulated 
with  a  high  fidelity  propagator  that  includes  nonlinear  effects,  Earth  oblateness,  and  navigation  and  control 
errors.  Burns  are  assumed  to  be  instantaneous  and  the  satellite  is  assumed  to  have  three  axis  burn 
capability.  Burn  calculations  are  corrupted  using  the  expected  distribution  of  the  thrust  magnitude,  which 
is  assumed  to  be  Gaussian. 

The  performance  of  the  satellite  will  be  evaluated  based  on  the  number  of  maneuvers  per  orbit,  the  amount 
of  AV  expended  per  day,  and  the  time  spent  within  the  deadband,  “availability”.  The  radius  of  the 
deadband  is  set  at  three  meters.  Using  this  deadband  and  running  short  simulations  without  Earth 
oblateness,  the  time  spent  during  a  maneuver  is  set  at  one  quarter  of  an  orbit  and  the  percentage  of  the 
deadband  at  which  the  maneuver  is  triggered  is  set  at  ninety  percent.  Using  this  percentage  of  the  deadband 
is  important  because  waiting  until  the  satellite  is  outside  of  the  deadband  before  a  maneuver  is  started  leads 
to  the  satellite  spending  more  of  its  time  outside  of  its  deadband  and  availability  suffers  greatly. 

For  the  Monte  Carlo  simulation,  the  following  parameters  in  Tables  1  and  2  are  used. 

Table  1 

BASELINE  ERROR  SOURCES 


Error  Source 

Position 

(mm) 

Velocity 

(mm/s) 

Thrust 

(%) 

Value  -  la 

5.0 

0.5 

5 
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Table  2 


STATIONKEEPING  PARAMETERS 


Constraint 

Deadband 

(m) 

Maneuver  Time 
(orbits) 

Threshold 

(%) 

Value 

3 

0.25 

90 

B.  In-plane  Case 

The  first  case  analyzed  is  the  in-plane  relative  formation  case.  For  initialization,  the  chief  satellite’s  initial 
mean  orbital  elements  will  be: 

a=  7378  km,  0=0,  i=50  deg,  qi=0,  q2=0,  Q=0  deg 
The  initial  mean  relative  orbital  elements  for  the  deputy  are: 

ae=  500  m,  xd=0,  yd=0,  zmax=0,  y=0  deg,  [HO  deg 

For  comparison  purposes,  a  baseline  case  is  run  using  only  the  unperturbed  two  body  equations  of  motion. 
The  results  from  this  simulation  are  the  optimal  results  achievable.  The  goal  of  the  guidance  approach  is  to 
get  the  formation  results  with  J2  perturbations  as  close  to  this  nominal  as  possible.  The  results  are  in  Table 
3  and  Figure  13  shows  a  histogram  of  the  number  of  maneuvers,  AV  per  day,  and  percent  of  time  spent 
within  the  deadband. 
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Figure  13:  Histogram  of  Baseline  Case 


Table  3 


BASELI1N 

IE  RESULTS 

Within 

Deadband 

AV 

(m/s/d  ay) 

Num.  Maneuvers 
per  orbit 

Mean 

99.0% 

0.264 

1.752 

Std.  Deviation 

0.9% 

0.077 

0.487 

To  gauge  the  efficiency  of  our  stationkeeping  approach,  the  Monte  Carlo  simulation  is  run  that  includes  J2 
in  the  simulation  dynamics,  but  not  in  the  guidance.  Instead,  it  employs  the  maneuver  scheme  used  in  a 
recent  AFRL  study15  (which  is  based  on  the  standard  unperturbed  dynamics).  As  can  be  seen  in  Table  4 
and  Figure  14,  the  results  from  this  simulation  are  heavily  affected  by  Earth  oblateness.  Both  the  amount 
of  AV  and  the  number  of  maneuvers  have  increased.  The  availability  is  also  significantly  reduced. 
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Figure  14:  Histogram  of  In-plane  Case  with  Earth  Oblateness 


Table  4 

RESULTS  WITH  EARTH  OBLATENESS 


Within 

Deadband 

AV 

(m/s/day) 

Num.  Maneuvers 
per  orbit 

Mean 

93.2% 

0.365 

2.308 

Std.  Deviation 

2.7% 

0.036 

0.152 

Finally,  for  the  in-plane  case,  we  include  our  maneuver  scheme  which  incorporates  the  J2  perturbations  into 
the  guidance  to  minimize  the  effect  of  Earth  oblateness.  The  previous  two  examples  for  the  in-plane  case 
used  simple  CWH  equations  in  their  maneuver  schemes  and  used  the  unperturbed  mean  motion.  It  can  be 
seen  from  the  results  in  Table  5  that  using  our  maneuver  scheme  and  using  the  perturbed  mean  motion 
allows  for  the  performance  of  this  scheme  to  be  greatly  improved  over  that  of  the  previous  case.  The 
performance  is  not  as  good  as  the  baseline,  but  the  guidance  approach  that  we  have  developed  greatly 
mitigates  the  impact  J2  has  on  the  formation  maintenance. 
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Figure  15:  Histogram  of  In-plane  Case  with  New  Guidance 


Table  5 

IN-PLANE  RESULTS  WITH  NEW  GUIDANCE 


Within 

Deadband 

AV 

(m/s/day) 

Num.  Maneuvers 
per  orbit 

Mean 

98.2% 

0.284 

1.829 

Std.  Deviation 

1 .4% 

0.085 

0.528 
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C.  Out-of-plane  Case 


For  the  out-of-plane  formation,  the  same  initial  conditions  are  used,  with  the  exception  that  zmax=500  for 
this  formation.  A  baseline  case  (with  the  simple  unperturbed  dynamic  model)  is  run  for  the  out-of-plane 
formation  and  the  results  are  in  Table  6  and  Figure  16.  The  results  are  very  similar  to  the  in-plane  case. 
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Figure  16:  Histogram  of  Baseline  Case 


Table  6 

OUT-OF-PLANE  BASELINE 


Within 

Deadband 

AV 

(m/s/d  ay) 

Num.  Maneuvers 
per  orbit 

Mean 

99.0% 

0.262 

1.746 

Std.  Deviation 

1.0% 

0.081 

0.512 
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The  next  case  run  uses  the  same  maneuver  scheme  as  the  baseline  but  now  we  include  J2  perturbations  in 
the  equations  of  motion  (but  not  in  the  guidance  scheme).  The  results  are  found  in  Table  7  and  Figure  17. 
The  performance  is  seriously  degraded  due  to  the  perturbations.  The  number  of  maneuvers  has 
significantly  increased  as  well  as  the  AV  and  the  availability  has  dropped  by  just  under  thirty  percent. 
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Figure  17:  Histogram  of  Out-of-Plane  with  Perturbations 


Table  7 

OUT-OF-PLANE  WITH  PERTURBATIONS 


Within 

Deadband 

AV 

(m/s/d  ay) 

Num.  Maneuvers 
per  orbit 

Mean 

71.9% 

0.537 

2.899 

Std.  Deviation 

4.3% 

0.038 

0.123 
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Next  the  new  guidance  scheme  with  J2  effects  is  employed  in  the  simulation.  The  results  are  improved  over 
the  previous  case.  Availability  is  increased  by  seventeen  percent  over  the  unmatched  case  and  the  AV 
requirements  and  the  number  of  maneuvers  are  decreased  by  thirty-five  percent  and  twenty-five  percent, 
respectively.  The  results  are  found  in  Table  8  and  Figure  18. 
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Figure  18:  Out-of-Plane  with  Oblateness  and  New  Guidance 


Table  8 

OUT-OF-PLANE  WITH  OBLATENESS  AND  NEW  GUIDANCE 


Within 

Deadband 

AV 

(m/s/day) 

Num.  Maneuvers 
per  orbit 

Mean 

89.2% 

0.348 

2.167 

Std.  Deviation 

3.5% 

0.073 

0.408 

To  assess  the  relative  merits  of  the  different  features  of  the  new  guidance  scheme,  we  run  a  case  without 
period  matching.  The  J2  effects  are  still  partially  accounted  for  in  this  case  by  using  the  perturbed  mean 
motion  in  the  calculation  of  the  nominal  trajectory.  The  results  are  shown  in  Table  9.  The  performance  is 
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better  than  the  case  with  no  J2  in  the  guidance  (Table  7)  but  worse  than  the  fully  developed  scheme  (Table 
8).  This  shows  that  both  period  matching  and  the  modifications  to  the  CWH  equations  are  important  in  the 
improvement  of  the  performance  over  previous  approaches. 


Table  9 


OUT-OF-PLANE  WITH  OBLATENESS  WITHOUT  PERIOD  MATCHING 


Within 

Deadband 

AV 

(m/s/day) 

Num.  Maneuvers 
per  orbit 

Mean 

87.9% 

0.379 

2.337 

Std.  Deviation 

4.3% 

0.064 

0.339 

D.  Summary  of  Results 

The  effects  of  J2  are  significant  for  relative  orbits  in  a  low  Earth  orbit.  Using  the  maneuver  scheme 
presented  in  this  report,  with  period  matching  and  using  a  perturbed  mean  motion  calculation,  brings  the 
performance  of  the  formation  back  close  to  the  baseline  case.  The  results  for  the  in-plane  case  with  the 
maneuver  scheme  are  very  close  to  that  of  the  baseline.  The  out-of-plane  results  are  improved  from  the 
case  with  perturbations  but  are  not  quite  at  baseline  due  to  the  drift  caused  by  the  differential  nodal  rate. 

The  process  presented  in  this  report  allows  for  very  tight  control  of  a  relative  trajectory  using  several 
stationkeeping  burns  per  orbit.  When  J2  is  accounted  for,  the  process  helps  to  minimize  the  amount  of  A  V 
required  to  keep  the  deputy  near  its  nominal  position. 
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CHAPTER  8.  FUTURE  WORK 


The  maneuver  scheme  laid  out  in  the  paper  helps  to  reduce  the  effects  of  J2  on  the  stationkeeping 
performance  of  a  satellite  maintaining  a  relative  trajectory.  However,  this  approach  has  not  quite  gotten  the 
number  of  maneuvers,  AV ,  and  availability  back  to  values  similar  to  the  baseline  cases.  The  next  step 
would  be  to  add  J2  guidance  into  state  transition  matrix  for  the  calculation  of  the  station-keeping  burn. 
Currently,  the  maneuver  state  transition  matrix  is  based  off  of  the  two  body  equations  of  motion  and  does 
not  take  into  account  J2.  Using  a  J2  relative  state  transition  matrix  should  help  to  move  the  simulation 
results  closer  to  the  baseline  cases.  An  attempt  was  made  to  accomplish  this  work,  but  the  current 
development  of  state  transition  matrices  with  J2  did  not  seem  to  have  the  sufficient  accuracy  for  this 
application.  This  area  merits  further  investigation. 

Another  next  step  in  the  research  would  be  to  account  for  uncertainty  and  eccentricity  in  the  chief 
trajectory.  Currently,  the  chief  is  presumed  to  be  perfectly  known  and  circular.  The  chief  will  have  some 
error  in  the  knowledge  of  its  position  and  velocity  and  this  will  also  affect  AV  and  the  maneuver  scheme. 
Eccentricity  in  the  chief  orbit  will  affect  the  relative  orbit  also. 
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APPENDIX  I.  PERIODIC  VARIATIONS  IN  ORBIT  ELEMENTS 


From  Ref.  10: 


0  =  (1-5  cos2  z)  1 


alp=  0 

e'p  =  - 

ilp  =  0 


^sin2  i'' 
8  a2 


(l- 100  cos2  i\ql  sin0  +  q2  cos  0) 


ip 

= 


ip 

q2  = 


f  qx  sin2  i  '' 
v  16a2  j 
(  _  -2  A 


(l- 100  cos2/) 


q2  sin  / 

76a2 


(l- 100  cos2/) 


J 


q!p  =  o 


aspl  =  ■ 


■  (l  -  3  cos2  / \q{  cos  6  +  q2  sin  0) 


2  aRt 

0spl  =  (l  -  5  cos2  / )(r/,  sin  #  -q2  cos  8) 

ispl  =  0 

sp  1 


Sol 

q2  = 


nsp'  = 


_  3(l — 3  cos  /)  ^2  CQS  Q  +  2qi+  qi  cos  2 0  +  q2  sin  2$) 
8a 

_  3(l  3  cos  0 ^ s jn  q  +  2 q2  +  ^  sin  26*  -  q2  cos  26) 
(g,  svt\0-q2  cos  6) 


8  a2 
9  cos  / 
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sp  2 

aF  =- 


3  sin2  i 


(l  +  3 qx  cos  0  +  3 q2  sin  6l)cos  20 


9sp  = 
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) 
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sin  29  -  Sm  y  (3  sin  9  -  sin  39) 
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flsp2  = - —  [3  sin  29  +  3 (ql  sin  9 +  q2  cos  9)  +  {q]  sin  29  -  q2  cos  36*)] 
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